#include "stdafx.h"
#include "FnSurfaceGlobalParam.h"
#include <cmath>
#include "./taucs/taucs_fix.h"
#include <set>
#include <queue>

int FnSurfaceGlobalParam::nx_face_side_idx = -1;
int FnSurfaceGlobalParam::px_face_side_idx = -1;
int FnSurfaceGlobalParam::ny_face_side_idx = -1;
int FnSurfaceGlobalParam::py_face_side_idx = -1;
int FnSurfaceGlobalParam::nz_face_side_idx = -1;
int FnSurfaceGlobalParam::pz_face_side_idx = -1;


FnSurfaceGlobalParam::FnSurfaceGlobalParam(Mesh& _mesh)
  : mesh_(_mesh), weight_type_(kUniform), side_idx_(-1)
{
  mesh().add_property(vprop_lns_idx_);

  mesh().add_property(vprop_param_);
  mesh().add_property(eprop_w_);
  mesh().add_property(fprop_side_idx_, "fprop_side_idx");
  mesh().add_property(vprop_side_idx_, "vprop_side_idx");
  mesh().add_property(vprop_type_, "vprop_type");
  mesh().add_property(eprop_type_, "eprop_type");  
  mesh().add_property(eprop_w_);
  mesh().add_property(eprop_edge_flag_);
  mesh().add_property(eprop_length_);
}

FnSurfaceGlobalParam::~FnSurfaceGlobalParam(void)
{
  mesh().remove_property(vprop_lns_idx_);

   mesh().remove_property(vprop_param_);
  mesh().remove_property(eprop_w_);
  mesh().remove_property(fprop_side_idx_);
  mesh().remove_property(vprop_side_idx_);
  mesh().remove_property(vprop_type_);
  mesh().remove_property(eprop_type_);  
  mesh().remove_property(eprop_w_);
  mesh().remove_property(eprop_length_);
  mesh().remove_property(eprop_edge_flag_);
	 
}


void FnSurfaceGlobalParam::update_edge_weights()
{
  Mesh::EdgeIter        e_it;
  Scalar                surface_w;

  for (e_it=mesh().edges_begin(); e_it!=mesh().edges_end(); ++e_it)
  {
  
    surface_w = 0;    
   // if (mesh().property(eprop_type_, e_it) == kSideEdge)
   // {
      if (weight_type_ == kUniform)         surface_w = 1;
      else if (weight_type_ == kCotangent)  surface_w = cot_surface_weight(e_it);
    //}
    
    mesh().property(eprop_w_, e_it) = surface_w;
  }
}

Scalar FnSurfaceGlobalParam::cot_triangle(const Mesh::HalfedgeHandle _corner) const
{
  Mesh::Point         p, q, r;
  OpenMesh::Vec3d  u, v;
  double              dot_val;

  p = mesh().point(mesh().to_vertex_handle(_corner));
  q = mesh().point(mesh().to_vertex_handle(mesh().next_halfedge_handle(_corner)));
  r = mesh().point(mesh().from_vertex_handle(_corner));

  u = OpenMesh::Vec3d(q[0]-p[0], q[1]-p[1], q[2]-p[2]);
  v = OpenMesh::Vec3d(r[0]-p[0], r[1]-p[1], r[2]-p[2]);

  dot_val = OpenMesh::dot(u, v);
  
  return  (Scalar) (dot_val / std::sqrt(u.sqrnorm()*v.sqrnorm() - dot_val*dot_val));
}



Scalar FnSurfaceGlobalParam::cot_surface_weight(const Mesh::EdgeHandle _eh) const
{
  if (mesh().property(eprop_type_, _eh) != kSideEdge)
    return  0;

  Mesh::FaceHandle                  fh;
  Scalar                            w;

  w = 0.5*(cot_triangle(mesh().halfedge_handle(_eh, 0))  + cot_triangle(mesh().halfedge_handle(_eh, 1)));

  return  w;
}


void FnSurfaceGlobalParam::init()
{
  Mesh::VertexIter            v_it;
  Mesh::EdgeIter              e_it;
  Mesh::FaceIter              f_it;
  Mesh::FaceVertexIter        fv_it;
  Mesh::HalfedgeHandle        heh[2];
  int                         lns_idx;

  Mesh::VertexEdgeIter  ve_it;  

  // set eprop_type_
  for (e_it=mesh().edges_begin(); e_it!=mesh().edges_end(); ++e_it)
  {
    mesh().property(eprop_type_, e_it) = kTypelessEdge;
    
    if(mesh().property(fprop_side_idx_, mesh().face_handle(mesh().halfedge_handle(e_it, 0)))
		== mesh().property(fprop_side_idx_, mesh().face_handle(mesh().halfedge_handle(e_it, 1)))){
      mesh().property(eprop_type_, e_it) = kSideEdge;
	}
    else 
      mesh().property(eprop_type_, e_it) = kWedgeEdge;
  }

  // set eprop_volume_w_ & eprop_surface_w_
  update_edge_weights();

  std::vector<Mesh::VertexHandle> corner_vertices;

  // set vprop_type_ & vprop_lns_idx_
  lns_idx = 0;
  for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
  {
    int         num_adj_edge_on_wedge = 0;
    int         num_adj_edge_on_side = 0;
    VertexType  vtype = kTypelessVertex;

    for (ve_it=mesh().ve_iter(v_it); ve_it; ++ve_it)
    {
	  if (mesh().property(eprop_type_, ve_it) == kSideEdge)         ++num_adj_edge_on_side;
      else if (mesh().property(eprop_type_, ve_it) == kWedgeEdge)   ++num_adj_edge_on_wedge;
    }

	mesh().property(vprop_side_idx_, v_it) = -1;
	if (num_adj_edge_on_wedge > 2){          
		vtype = kCornerVertex;
	}
    else if (num_adj_edge_on_wedge == 2)    vtype = kWedgeVertex;
	else if (num_adj_edge_on_side > 0){      
		vtype = kSideVertex;
		mesh().property(vprop_side_idx_, v_it) = mesh().property(fprop_side_idx_, mesh().face_handle(mesh().halfedge_handle(v_it)));
	}

    mesh().property(vprop_type_, v_it) = vtype;
    mesh().property(vprop_lns_idx_, v_it) = (vtype != kCornerVertex) ? lns_idx++ : -1;

	if(vtype == kCornerVertex){
		corner_vertices.push_back(v_it);
	}
  }

  corner_infos_.clear();
  corner_infos_.resize(corner_vertices.size());

  Mesh::VertexFaceIter vf_it;
  int i;

  for(i = 0; i < corner_infos_.size(); i++){
	  corner_infos_[i].corner_vtx_ = corner_vertices[i];
	  for(vf_it = mesh().vf_iter(corner_infos_[i].corner_vtx_); vf_it; ++vf_it){
		  corner_infos_[i].nbhd_patch_idx_.insert(mesh().property(fprop_side_idx_, vf_it));
	  }

	  assert(corner_infos_[i].nbhd_patch_idx_.size() >= 3);
  }

 /* // set linear system index (i.e., vprop_lns_idx_)
  for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
  {
    mesh().property(vprop_lns_idx_, v_it) = -1;
  }
  
  lns_idx = 0;
  for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
  {
    if (mesh().property(vprop_type_, v_it) != kCornerVertex)    
      mesh().property(vprop_lns_idx_, v_it) = lns_idx++;
  }
*/
   // give a patch idx to each vertex
  // and set patch infos (corners & adj patches)
  set_patch_infos();

}

void FnSurfaceGlobalParam::update_patch_paths()
{
	int i, k, j;
	std::vector<Mesh::HalfedgeHandle> path;
	Mesh::EdgeIter e_it;
	Mesh::FaceIter f_it;
	Mesh::FaceEdgeIter fe_it;
	Mesh::Point           pnt[2];
	float max_edge_length_ = 0.0f, len;
	OpenMesh::FPropHandleT<bool> fprop_flag;

	std::vector<Mesh::FaceHandle> faces;
	std::vector<Mesh::EdgeHandle> edges;

	OpenMesh::EPropHandleT<bool> eprop_cube_edge_flag;

	mesh().get_property_handle(eprop_cube_edge_flag, "eprop_cube_edge_flag");

	
	mesh().add_property(fprop_flag);
	
	 for(e_it=mesh().edges_begin(); e_it != mesh().edges_end(); ++e_it){
	    pnt[0] = mesh().point( mesh().to_vertex_handle(mesh().halfedge_handle(e_it, 0)));
        pnt[1] = mesh().point( mesh().to_vertex_handle(mesh().halfedge_handle(e_it, 1)));

		len = (pnt[0]-pnt[1]).length();

		mesh().property(eprop_length_, e_it) = len;
		mesh().property(eprop_edge_flag_, e_it) = true;
		mesh().property(eprop_cube_edge_flag, e_it) = true;

		max_edge_length_ += len;
	 }

	 for(f_it = mesh().faces_begin(); f_it != mesh().faces_end(); ++f_it){
		 if(mesh().property(fprop_side_idx_, f_it) == kNone || mesh().property(fprop_side_idx_, f_it) == kCorner){
			 if(mesh().property(fprop_side_idx_, f_it) == kCorner)
				 mesh().property(fprop_side_idx_, f_it) = kNone;

			 for(fe_it = mesh().fe_iter(f_it); fe_it; ++fe_it){
				 mesh().property(eprop_edge_flag_, fe_it) = false;
			 }
		 }

		 mesh().property(fprop_flag, f_it) = false;
	 }

	 Mesh::FaceHalfedgeIter fh_it;
	 Mesh::FaceHandle fh, fh2;
	 

	for(i = 0; i < patch_infos_.size(); i++){
		std::queue<Mesh::FaceHandle> queue;
		faces.clear();
		edges.clear();
	 
		for(k = 0; k < 4; k++){
			get_patch_path(patch_infos_[i].corners_[k], patch_infos_[i].corners_[(k+1)%4], path);

			for(j = 0; j < path.size(); j++){
				mesh().property(eprop_edge_flag_, mesh().edge_handle(path[j])) = true;
				edges.push_back(mesh().edge_handle(path[j]));
				fh2 = mesh().face_handle(path[j]);

				fh = mesh().face_handle(mesh().opposite_halfedge_handle(path[j]));
				mesh().property(fprop_flag, fh) = true;
				faces.push_back(fh);
				
				if(mesh().property(fprop_side_idx_, fh2) != kNone)
					continue;

				mesh().property(fprop_side_idx_, fh2) = i;

				//mesh().property(eprop_edge_flag_, mesh().edge_handle(path[j])) = false;
		
				for(fh_it = mesh().fh_iter(fh2); fh_it; ++fh_it){
					if(Mesh::HalfedgeHandle(fh_it) == path[j]){
						continue;
					}

					fh = mesh().face_handle(mesh().opposite_halfedge_handle(fh_it));

					if(mesh().property(fprop_side_idx_, fh) == kNone && mesh().property(fprop_flag, fh) == false){
						queue.push(fh);
					}
				}
			}

		
		}
		
		while(!queue.empty()){
			fh = queue.front();
			queue.pop();
			
			if(mesh().property(fprop_side_idx_, fh) == kNone && mesh().property(fprop_flag, fh) == false){
				mesh().property(fprop_side_idx_, fh) = i;
				for(fh_it = mesh().fh_iter(fh); fh_it; ++fh_it){
				
					fh2 = mesh().face_handle(mesh().opposite_halfedge_handle(fh_it));
					if(mesh().property(fprop_side_idx_, fh2) == kNone  && mesh().property(fprop_flag, fh) == false){
			//			mesh().property(fprop_side_idx_, fh2) = i;
						queue.push(fh2);
					}
				}
			}
		}


		for(k = 0; k < faces.size(); k++){
			mesh().property(fprop_flag, faces[k]) = false;
		}

		for(k = 0; k < edges.size(); k++){
			mesh().property(eprop_edge_flag_, edges[k]) = false;
			mesh().property(eprop_cube_edge_flag, edges[k]) = false;
		}
	}

	mesh().remove_property(fprop_flag);

	int idx_f, idx_e;
	int cnt;
	std::set<int> indices;
	for(f_it = mesh().faces_begin(); f_it != mesh().faces_end(); ++f_it){
		idx_f = mesh().property(fprop_side_idx_, f_it);
		cnt = 0;
		indices.clear();

		for(fh_it = mesh().fh_iter(f_it); fh_it; ++fh_it){
			fh = mesh().face_handle(mesh().opposite_halfedge_handle(fh_it));
			idx_e = mesh().property(fprop_side_idx_, fh);

			if(idx_e != idx_f){
				cnt++;
				indices.insert(idx_e);
			}
		}

		if(indices.size() == 1 && cnt >= 2){
			mesh().property(fprop_side_idx_, f_it) = idx_e;
		}
	}
}

bool FnSurfaceGlobalParam::get_patch_path(const Mesh::VertexHandle _start_vh, const Mesh::VertexHandle _end_vh, std::vector<Mesh::HalfedgeHandle> &path)
{
	OpenMesh::VPropHandleT<Mesh::HalfedgeHandle> VPropHHandle;
	OpenMesh::VPropHandleT<float> VPropPriority;
	OpenMesh::VPropHandleT<bool> VPropAuxFlag;
	OpenMesh::VPropHandleT<bool> VPropInQueue;

	
	 mesh().add_property(VPropPriority);
	 mesh().add_property(VPropHHandle);
	 mesh().add_property(VPropAuxFlag);
	 mesh().add_property(VPropInQueue);

	 Mesh::EdgeIter e_it;
	 Mesh::VertexIter	  v_it;
     
	 std::priority_queue<VertexPriority> pq;
	 float len;

	 for(v_it=mesh().vertices_begin(); v_it != mesh().vertices_end(); ++v_it){
		 mesh().property(VPropPriority, v_it) = max_edge_length_;
		 mesh().property(VPropAuxFlag, v_it) = false;
		 mesh().property(VPropInQueue, v_it) = false;
	 }

	 mesh().property(VPropPriority, _start_vh) = 0.0f;

	 VertexPriority pr;
	 pr.priority = 0.0f;
	 pr.vh_ = _start_vh;

	 Mesh::VertexHandle _vh, _nvh;

	 mesh().property(VPropInQueue, _start_vh) = true;

	   //CString str;
		

	 pq.push(pr);

	   //str.Format(_T("push [%d] : %f\n"), _start_vh.idx(), property(VPropPriority, _start_vh));

		//OutputDebugString(str);

	 while(!pq.empty()){

		do{
	       pr = pq.top();
		   pq.pop();
		   _vh = pr.vh_;

		  //str.Format(_T("pop [%d] : %f\n"), _vh.idx(), property(VPropPriority, _vh));

		//OutputDebugString(str);

		}while(mesh().property(VPropInQueue, _vh) == false && !pq.empty());

		/*if(pq.empty())
			break;*/

		if(mesh().property(VPropInQueue, _vh) == false){
			break;
		}

		mesh().property(VPropInQueue, _vh) = false;

		mesh().property(VPropAuxFlag, _vh) = true;

		if(_vh == _end_vh)
			break;


			Mesh::VertexOHalfedgeIter vohe_it;

			for(vohe_it=mesh().voh_iter(_vh); vohe_it; ++vohe_it){
				_nvh = mesh().to_vertex_handle(vohe_it);

				len = mesh().property(eprop_length_, mesh().edge_handle(vohe_it));

				if(mesh().property(VPropAuxFlag, _nvh))
					continue;

				if(mesh().property(eprop_edge_flag_, mesh().edge_handle(vohe_it)) == true)
					continue;

				if(mesh().property(fprop_side_idx_, mesh().face_handle(vohe_it)) != kNone && mesh().property(fprop_side_idx_, mesh().face_handle(mesh().opposite_halfedge_handle(vohe_it))) != kNone
					&& mesh().property(fprop_side_idx_, mesh().face_handle(vohe_it)) == mesh().property(fprop_side_idx_, mesh().face_handle(mesh().opposite_halfedge_handle(vohe_it))))
					continue;

				if(mesh().property(VPropInQueue, _nvh) && mesh().property(VPropPriority, _nvh) > mesh().property(VPropPriority, _vh) + len){
					mesh().property(VPropInQueue, _nvh) = false;
				}

				if( !mesh().property(VPropInQueue, _nvh) ){
					mesh().property(VPropPriority, _nvh) = mesh().property(VPropPriority, _vh) + len;
					mesh().property(VPropHHandle, _nvh) = vohe_it;

					mesh().property(VPropInQueue, _nvh) = true;

					pr.priority = mesh().property(VPropPriority, _vh) + len;
					pr.vh_ = _nvh;
					pq.push(pr);

					//str.Format(_T("push [%d] : %f\n"), _nvh.idx(), property(VPropPriority, _nvh));

					//OutputDebugString(str);
				}
			}
	 }

	 path.clear();

	_vh = _end_vh;
	Mesh::HalfedgeHandle _he;
	bool is_path_valid = true;

	do{
		_he = mesh().property(VPropHHandle, _vh);

		if(_he.is_valid() == false){
			is_path_valid = false;
			path.clear();
			break;
		}

		_vh = mesh().from_vertex_handle(_he);
		path.insert(path.begin(), _he);

	}while(_vh != _start_vh);

 	 mesh().remove_property(VPropPriority);
	 mesh().remove_property(VPropHHandle);
	 mesh().remove_property(VPropAuxFlag);
	 mesh().remove_property(VPropInQueue);

	return is_path_valid;
}

void FnSurfaceGlobalParam::update_patches()
{
	 swap_param_values();
 
	 Mesh::FaceIter f_it;
	 Mesh::FaceVertexIter fv_it;
	 std::set<int> indices;
	 int idx, i, j;
	 Mesh::FaceHalfedgeIter fh_it;
	 TransFunc tf;
	 Param new_uv[3], uv;

	 std::vector<Mesh::VertexHandle> new_corners, old_corners;
	 std::vector<int> corner_indices;
	 Mesh::VertexVertexIter vv_it;
	 
	
	 for(f_it = mesh().faces_begin(); f_it != mesh().faces_end(); ++f_it){

		 indices.clear();
		 bool corner = false;


		 for(fv_it = mesh().fv_iter(f_it); fv_it; ++fv_it){
			 if(mesh().property(vprop_type_, fv_it) == kCornerVertex){
				 corner = true;
				 continue;
			 }

			 idx = mesh().property(vprop_side_idx_, fv_it);
			 indices.insert(idx);
		 }

		 if(corner == false && indices.size() == 1){
			 mesh().property(fprop_side_idx_, f_it) = idx;
		 }
		 else{
			 mesh().property(fprop_side_idx_, f_it) = kNone;

			 if(corner == false && indices.size() == 3){
				 mesh().property(fprop_side_idx_, f_it) = kCorner;

				 std::vector<int> indices2;

				 int idx;
				 Mesh::VertexHandle corner_vh;

				 for(i = 0; i < corner_infos_.size(); i++){
					 //if(indices == corner_infos_[i].nbhd_patch_idx_){
					 if(std::includes(corner_infos_[i].nbhd_patch_idx_.begin(), corner_infos_[i].nbhd_patch_idx_.end(),
						 indices.begin(), indices.end())){
						 idx = i;
						 //corner_indices.push_back(i);
						 corner_vh = corner_infos_[i].corner_vtx_;
						 //old_corners.push_back(corner_vh);

						//tuan
						 mesh().property(vprop_type_, corner_vh) = kTypelessVertex;
						 for(vv_it = mesh().vv_iter(corner_vh); vv_it; ++vv_it){
							 if(mesh().property(vprop_side_idx_, vv_it) >= 0){
								 mesh().property(vprop_side_idx_, corner_vh) = mesh().property(vprop_side_idx_, vv_it);
								 break;
							 }
						 }
						  //

						 indices2.resize(corner_infos_[i].nbhd_patch_idx_.size());
						 std::copy(corner_infos_[i].nbhd_patch_idx_.begin(), corner_infos_[i].nbhd_patch_idx_.end(), 
							 indices2.begin());

						 break;
					 }
				 }

				 fv_it = mesh().fv_iter(f_it);
				 //new_corners.push_back(fv_it);

				 //tuan
				 mesh().property(vprop_type_,fv_it) = kCornerVertex;
				 mesh().property(vprop_side_idx_, fv_it) = -1;
				 corner_infos_[idx].corner_vtx_ = fv_it;
				 //

				 for(i = 0; i < indices2.size(); i++){
					 for(j = 0; j < 4; j++){
						 if(patch_infos_[indices2[i]].corners_[j] == corner_vh){
							 patch_infos_[indices2[i]].corners_[j] = fv_it;
						 }
					 }
				 }

			 }	

		 }
	 }

	 if (new_corners.size() > nSides + 2)
	 {
		 printf("number of corners is > %d\n", nSides + 2);
	 }

	 for(int i = 0; i < new_corners.size(); i++){
		 mesh().property(vprop_type_, new_corners[i]) = kCornerVertex;
		 mesh().property(vprop_side_idx_, new_corners[i]) = -1;
	 }
	 
	 for(int i = 0; i < old_corners.size(); i++){
		 mesh().property(vprop_type_, old_corners[i]) = kTypelessVertex;

		 for(vv_it = mesh().vv_iter(old_corners[i]); vv_it; ++vv_it){
			 if(mesh().property(vprop_side_idx_, vv_it) >= 0){
				 mesh().property(vprop_side_idx_, old_corners[i]) = mesh().property(vprop_side_idx_, vv_it);
				 break;
			 }
		 }
	 }
	 for(int j = 0; j < corner_indices.size(); j++){
		 corner_infos_[corner_indices[j]].corner_vtx_ = new_corners[j];
	 }
	 // mesh().remove_property(fprop_flags);

	 for(f_it = mesh().faces_begin(); f_it != mesh().faces_end(); ++f_it){

		 indices.clear();
		 bool corner = false;


		 for(fv_it = mesh().fv_iter(f_it); fv_it; ++fv_it){
			 if(mesh().property(vprop_type_, fv_it) == kCornerVertex){
				 corner = true;
				 continue;
			 }

			 idx = mesh().property(vprop_side_idx_, fv_it);
			 indices.insert(idx);
		 }

		 if(corner == false && indices.size() == 1){
			 mesh().property(fprop_side_idx_, f_it) = idx;
		 }

	 }


	 int lns_idx;

	 Mesh::VertexIter v_it;
	 for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
	 {
		 mesh().property(vprop_lns_idx_, v_it) = -1;
	 }

	 lns_idx = 0;
	 for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
	 {
		 if (mesh().property(vprop_type_, v_it) != kCornerVertex)    
			 mesh().property(vprop_lns_idx_, v_it) = lns_idx++;
	 }

}

void FnSurfaceGlobalParam::calc()
{
  //export_linear_system("a.dat", "b.dat");
  calc_with_taucs();
  //update_patches();
}

void FnSurfaceGlobalParam::save(const std::string& _filename)
{
  FILE    *fp;
  int     vidx;

  fp = fopen(_filename.c_str(), "w");

  Mesh::VertexIter           v_it;
  Mesh::Point           p;
  Param                 uv;
  OpenMesh::Vec2f    st;

  OpenMesh::VPropHandleT<int>  vprop_new_idx, vprop_uv_idx;
  OpenMesh::VPropHandleT<bool> vprop_flags;

  mesh().add_property(vprop_new_idx);
  mesh().add_property(vprop_uv_idx);
  mesh().add_property(vprop_flags);

  
  // v
  vidx = 0;
  for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
  {
  
    mesh().property(vprop_new_idx, v_it) = vidx++;
    p = mesh().point(v_it);

	fprintf(fp, "v %f %f %f\n", p[0], p[1], p[2]);
  }

  fprintf(fp, "vt 0 0\n");
  fprintf(fp, "vt 1 0\n");
  fprintf(fp, "vt 1 1\n");
  fprintf(fp, "vt 0 1\n");
  // vt
  vidx = 4;

  /*for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
  {
	if(mesh().property(vprop_lns_idx_, v_it) < 0)
		  continue;

	mesh().property(vprop_uv_idx, v_it) = vidx++;
    uv = mesh().property(vprop_param_, v_it);

    fprintf(fp, "vt %f %f\n", uv[0], uv[1]);
  }*/

  Mesh::FaceIter f_it;
  Mesh::FaceVertexIter fv_it;
  Mesh::FaceHandle fh;

  std::vector<std::vector<Mesh::FaceHandle> > patch_faces;
  std::vector<Mesh::VertexHandle> vertices;

  patch_faces.resize(6);

  int p_idx;
  for(f_it=mesh().faces_begin(); f_it!=mesh().faces_end(); ++f_it){
	p_idx = mesh().property(fprop_side_idx_, f_it);
	patch_faces[p_idx].push_back(f_it);
  }

  for(v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it){
	  mesh().property(vprop_flags, v_it) = false;
  }


  for(int i = 0; i < 6; i++){
	  vertices.clear();

	  for(int k = 0; k < patch_faces[i].size(); k++){
		fh = patch_faces[i][k];

		for(fv_it=mesh().fv_iter(fh); fv_it; ++fv_it){
			if(mesh().property(vprop_lns_idx_, fv_it) < 0)
				  continue;

			if(mesh().property(vprop_flags, fv_it) == false){
				vertices.push_back(fv_it);
				mesh().property(vprop_flags, fv_it) = true;
			}
		}

	  }

	  for(int k = 0; k < vertices.size(); k++){
		  mesh().property(vprop_flags, vertices[k]) = false;

		  if(mesh().property(vprop_side_idx_, vertices[k]) != i){
			  //assert(mesh().property(vprop_type_, vertices[k]) == kWedgeVertex);

			  swap_param_value(vertices[k]);
		  }
		  mesh().property(vprop_uv_idx, vertices[k]) = vidx++;
		  uv = mesh().property(vprop_param_, vertices[k]);
		  fprintf(fp, "vt %f %f\n", uv[0], uv[1]);
	  }

	  
	 mesh().property(vprop_uv_idx, patch_infos_[i].corners_[0]) = 0;
	 mesh().property(vprop_uv_idx, patch_infos_[i].corners_[1]) = 1;
	 mesh().property(vprop_uv_idx, patch_infos_[i].corners_[2]) = 2;
	 mesh().property(vprop_uv_idx, patch_infos_[i].corners_[3]) = 3;

	 
	 for(int k = 0; k < patch_faces[i].size(); k++){
	 fh = patch_faces[i][k];

	  fprintf(fp, "f ");
	  for(fv_it=mesh().fv_iter(fh); fv_it; ++fv_it){
		  fprintf(fp, "%d/%d ", mesh().property(vprop_new_idx, fv_it)+1, mesh().property(vprop_uv_idx, fv_it)+1);
	  }
	  fprintf(fp, "\n");
	 }
  }
 
  mesh().remove_property(vprop_new_idx);
  mesh().remove_property(vprop_uv_idx);
  mesh().remove_property(vprop_flags);

}





void FnSurfaceGlobalParam::export_linear_system(const std::string& _A_dat, const std::string& _b_dat)
{
	Mesh::VertexIter v_it;
	Mesh::VertexOHalfedgeIter voh_it;
	Mesh::VertexHandle vh_i, vh_j;

	FILE *fp_A, *fp_B;

  fp_A = fopen("a.dat", "w");
  fp_B = fopen("b.dat", "w");

	  // give a patch idx to each vertex
  // and set patch infos (corners & adj patches)
 set_patch_infos();

 std::vector<std::vector<Mesh::VertexHandle> > patch_vertices;
 int p_idx, p, i, j, k;

 patch_vertices.resize(6);

 for(v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it){
	 if(mesh().property(vprop_lns_idx_, v_it) < 0)
		 continue;

	 p_idx = mesh().property(vprop_side_idx_, v_it);
	 patch_vertices[p_idx].push_back(v_it);
 }

 Scalar w_ii, w_ij, bu, bv;
 int adj_idx;

 for(p = 0; p < 6; p++){
	 mesh().property(vprop_param_, patch_infos_[p].corners_[0]) = Param(0, 0);
	 mesh().property(vprop_param_, patch_infos_[p].corners_[1]) = Param(1, 0);
	 mesh().property(vprop_param_, patch_infos_[p].corners_[2]) = Param(1, 1);
	 mesh().property(vprop_param_, patch_infos_[p].corners_[3]) = Param(0, 1);

	 for(k = 0; k < patch_vertices[p].size(); k++){
		vh_i = patch_vertices[p][k];
		i = mesh().property(vprop_lns_idx_, vh_i);

		w_ii = 0.0;
		bu = 0.0;
		bv = 0.0;
		
		for(voh_it = mesh().voh_iter(vh_i); voh_it; ++voh_it){
			vh_j = mesh().to_vertex_handle(voh_it);
			j = mesh().property(vprop_lns_idx_, vh_j);

			w_ij = mesh().property(eprop_w_, mesh().edge_handle(voh_it));
			w_ii += w_ij;

			if(j >= 0){
				assert(mesh().property(vprop_side_idx_, vh_i) == p);
				adj_idx = mesh().property(vprop_side_idx_, vh_j);

				if(p != adj_idx){
					if(patch_infos_[p].adj_patch_idx_[0] == adj_idx){ // gamma -> alpha
						if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[3]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
							bv -= w_ij;
						}
						else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[2]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, -w_ij);
							bu += w_ij;
							bv -= w_ij;
						}
						else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[1]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, w_ij);
							bu += w_ij;
						}
						else{
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, w_ij);
						}
					}
					else if(patch_infos_[p].adj_patch_idx_[1] == adj_idx){ // beta -> alpha
						if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[0]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
							bu += w_ij;
						}
						else if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[3]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, -w_ij);
							bu += 2*w_ij;
						}
						else if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[2]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, w_ij);
							bu += 2*w_ij;
							bv += w_ij;
						}
						else{
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, w_ij);
							bu += w_ij;
							bv += w_ij;
						}
					}
					else if(patch_infos_[p].adj_patch_idx_[2] == adj_idx){ // alpha -> gamma
						if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[0]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
							bv += w_ij;
						}
						else if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[3]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, -w_ij);
							bu += w_ij;
							bv += w_ij;
						}
						else if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[2]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, w_ij);
							bu += w_ij;
							bv += 2*w_ij;
						}
						else{
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, w_ij);
							bv += 2*w_ij;
						}
					}
					else{ // alpha -> beta
						assert(patch_infos_[p].adj_patch_idx_[3] == adj_idx);
						if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[1]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
							bu -= w_ij;
						}
						else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[0]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, -w_ij);
						}
						else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[3]){
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, w_ij);
							bv += w_ij;
						}
						else{
							fprintf(fp_A, "%d %d %f\n", 1+2*i, 2+2*j, -w_ij);
							fprintf(fp_A, "%d %d %f\n", 2+2*i, 1+2*j, w_ij);
							bu -= w_ij;
							bv += w_ij;
						}
					}
				}
				else{
					fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
					fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
				}
			}
			else{
				Param uv = mesh().property(vprop_param_, vh_j);

				bu += w_ij*uv[0];
				bv += w_ij*uv[1];

				assert(vh_j == patch_infos_[p].corners_[0] || vh_j == patch_infos_[p].corners_[1] ||
					vh_j == patch_infos_[p].corners_[2] || vh_j == patch_infos_[p].corners_[3]);
			}
		}
		fprintf(fp_A, "%d %d %f\n", 1+2*i, 1+2*i, w_ii);
		fprintf(fp_A, "%d %d %f\n", 2+2*i, 2+2*i, w_ii);
		fprintf(fp_B, "%d %f\n", 1+2*i, bu);
		fprintf(fp_B, "%d %f\n", 2+2*i, bv);
	 }
	 
 }

  fclose(fp_A);
  fclose(fp_B);
}



void FnSurfaceGlobalParam::set_patch_infos()
{

	Mesh::VertexIter v_it;
	Mesh::VertexOHalfedgeIter  voh_it;
	int i, j;

	std::vector<Mesh::VertexHandle> corners;

	patch_infos_.clear();
	patch_infos_.resize(nSides);

	OpenMesh::VPropHandleT<bool> vflags;
	mesh().add_property(vflags);


	for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)  
	{
		mesh().property(vflags, v_it) = false;

		i = mesh().property(vprop_lns_idx_, v_it);
		if (i < 0){
			corners.push_back(v_it);
		}
		else if(mesh().property(vprop_type_, v_it) == kSideVertex){
			int idx = mesh().property(fprop_side_idx_, mesh().face_handle(mesh().halfedge_handle(v_it)));

			mesh().property(vprop_side_idx_, v_it) = idx;

		}
	}

  Mesh::VertexHandle vh;
  for(i = 0; i < corners.size(); i++){
	  vh = corners[i];

	  for(Mesh::VertexFaceIter vf_it = mesh().vf_iter(vh); vf_it; ++vf_it){
		  int side_itx = mesh().property(fprop_side_idx_, vf_it);
		
		  patch_infos_[side_itx].corners_[0] = vh;
	  }
  }

  for(i = 0; i < nSides; i++){
	  for(j = 0; j < 4; j++){
		  patch_infos_[i].adj_patch_idx_[j] = -1;
	  }
  }

  Mesh::VertexHandle corner_vh, vh_i;
  
  
  int pn, cnt;
  for(pn = 0; pn < patch_infos_.size(); pn++){
	  vh = corner_vh = patch_infos_[pn].corners_[0];
	  cnt = 1;
	  
	  do{
		  for(voh_it = mesh().voh_iter(vh); voh_it; ++voh_it){
			  if(mesh().property(eprop_type_, mesh().edge_handle(voh_it)) == kSideEdge)
				  continue;

			  if(mesh().property(fprop_side_idx_, mesh().face_handle(voh_it)) == pn){
				  vh = mesh().to_vertex_handle(voh_it);

				  if(patch_infos_[pn].adj_patch_idx_[cnt-1] == -1){
					  patch_infos_[pn].adj_patch_idx_[cnt-1] = mesh().property(fprop_side_idx_, mesh().face_handle(mesh().opposite_halfedge_handle(voh_it)));
				  }

				  if(mesh().property(vflags, vh) == false && mesh().property(vprop_lns_idx_, vh) >= 0){
					  mesh().property(vflags, vh) = true;
					  mesh().property(vprop_side_idx_, vh) = pn;
				  }
				  if(mesh().property(vprop_lns_idx_, vh) < 0 && vh != corner_vh){ // corner vertex
					  patch_infos_[pn].corners_[cnt] = vh;
				  }
				  break;
			  }
		  }

		  if(mesh().property(vprop_type_, vh) == kCornerVertex && vh != corner_vh){
			  cnt++;
		  }
	  }while(vh != corner_vh);

  }

  mesh().remove_property(vflags);

  int adj_idx;
  for(int p = 0; p < nSides; p++){
	  adj_idx = patch_infos_[p].adj_patch_idx_[0];// gamma -> alpha

	  if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[3]){
		  patch_infos_[p].trans_func_[0].a = 1.0;
		  patch_infos_[p].trans_func_[0].b = 0.0;
		  patch_infos_[p].trans_func_[0].t = Param(0.0, 1.0);
	  }
	  else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[2]){
		  patch_infos_[p].trans_func_[0].a = 0.0;
		  patch_infos_[p].trans_func_[0].b = 1.0;
		  patch_infos_[p].trans_func_[0].t = Param(1.0, 1.0);
	  }
	  else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[1]){
		  patch_infos_[p].trans_func_[0].a = -1.0;
		  patch_infos_[p].trans_func_[0].b = 0.0;
		  patch_infos_[p].trans_func_[0].t = Param(1.0, 0.0);
	  }
	  else{
		  patch_infos_[p].trans_func_[0].a = 0.0;
		  patch_infos_[p].trans_func_[0].b = -1.0;
		  patch_infos_[p].trans_func_[0].t = Param(0.0, 0.0);
	  }

	  adj_idx = patch_infos_[p].adj_patch_idx_[1];// beta -> alpha

	  if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[0]){
		  patch_infos_[p].trans_func_[1].a = 1.0;
		  patch_infos_[p].trans_func_[1].b = 0.0;
		  patch_infos_[p].trans_func_[1].t = Param(-1.0, 0.0);
	  }
	  else if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[3]){
		  patch_infos_[p].trans_func_[1].a = 0.0;
		  patch_infos_[p].trans_func_[1].b = 1.0;
		  patch_infos_[p].trans_func_[1].t = Param(0.0, 2.0);
	  }
	  else if(patch_infos_[p].corners_[1] == patch_infos_[adj_idx].corners_[2]){
		  patch_infos_[p].trans_func_[1].a = -1.0;
		  patch_infos_[p].trans_func_[1].b = 0.0;
		  patch_infos_[p].trans_func_[1].t = Param(2.0, 1.0);
	  }
	  else{
		  patch_infos_[p].trans_func_[1].a = 0.0;
		  patch_infos_[p].trans_func_[1].b = -1.0;
		  patch_infos_[p].trans_func_[1].t = Param(1.0, -1.0);
	  }

	  adj_idx = patch_infos_[p].adj_patch_idx_[2]; // alpha -> gamma
	  if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[0]){
		  patch_infos_[p].trans_func_[2].a = 1.0;
		  patch_infos_[p].trans_func_[2].b = 0.0;
		  patch_infos_[p].trans_func_[2].t = Param(0.0, -1.0);
	  }
	  else if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[3]){
		  patch_infos_[p].trans_func_[2].a = 0.0;
		  patch_infos_[p].trans_func_[2].b = 1.0;
		  patch_infos_[p].trans_func_[2].t = Param(-1.0, 1.0);
	  }
	  else if(patch_infos_[p].corners_[3] == patch_infos_[adj_idx].corners_[2]){
		  patch_infos_[p].trans_func_[2].a = -1.0;
		  patch_infos_[p].trans_func_[2].b = 0.0;
		  patch_infos_[p].trans_func_[2].t = Param(1.0, 2.0);
	  }
	  else{
		  patch_infos_[p].trans_func_[2].a = 0.0;
		  patch_infos_[p].trans_func_[2].b = -1.0;
		  patch_infos_[p].trans_func_[2].t = Param(2.0, 0.0);
	  }

	  adj_idx = patch_infos_[p].adj_patch_idx_[3]; // alpha -> beta

	  if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[1]){
		  patch_infos_[p].trans_func_[3].a = 1.0;
		  patch_infos_[p].trans_func_[3].b = 0.0;
		  patch_infos_[p].trans_func_[3].t = Param(1.0, 0.0);
	  }
	  else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[0]){
		  patch_infos_[p].trans_func_[3].a = 0.0;
		  patch_infos_[p].trans_func_[3].b = 1.0;
		  patch_infos_[p].trans_func_[3].t = Param(0.0, 0.0);
	  }
	  else if(patch_infos_[p].corners_[0] == patch_infos_[adj_idx].corners_[3]){
		  patch_infos_[p].trans_func_[3].a = -1.0;
		  patch_infos_[p].trans_func_[3].b = 0.0;
		  patch_infos_[p].trans_func_[3].t = Param(0.0, 1.0);
	  }
	  else{
		  patch_infos_[p].trans_func_[3].a = 0.0;
		  patch_infos_[p].trans_func_[3].b = -1.0;
		  patch_infos_[p].trans_func_[3].t = Param(1.0, 1.0);
	  }

  }
  printf("extract patch info ends...\n");
}

int FnSurfaceGlobalParam::calc_with_taucs()
{
	taucs_ccs_matrix  *A = NULL;
	taucs_double      *u;
	taucs_double      *bu;
	FILE*             fp_ijv;
	int               i, j, p, k;
	int               n;
	int               rc;
	char*             ijv_filename = "taucs_matrix.ijv";

	printf("FnSurfaceHarmonicParam::calc_with_taucs()\n");

	Mesh::VertexIter            v_it;
	Mesh::VertexHandle          vh_i, vh_j;
	Mesh::VertexOHalfedgeIter	voh_it;

	// calc # of free vertices (i.e., n x n)
	n = 0;
	for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
	{
		if (mesh().property(vprop_lns_idx_, v_it) < 0)
			continue;
		++n;
	}
	printf("n=%d\n", n);


	// fill the matrix A and bu, bv;
	fp_ijv = fopen(ijv_filename, "w");
	bu = (taucs_double *) taucs_vec_create(n*2, TAUCS_DOUBLE);


	std::vector<std::vector<Mesh::VertexHandle> > patch_vertices;
	int p_idx;

	patch_vertices.resize(nSides);

	for(v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it){
		if(mesh().property(vprop_lns_idx_, v_it) < 0)
			continue;

		p_idx = mesh().property(vprop_side_idx_, v_it);
		patch_vertices[p_idx].push_back(v_it);
	}

	Scalar w_ii, w_ij;
	int adj_idx;

	for(p = 0; p < nSides; p++){
		mesh().property(vprop_param_, patch_infos_[p].corners_[0]) = Param(0, 0);
		mesh().property(vprop_param_, patch_infos_[p].corners_[1]) = Param(1, 0);
		mesh().property(vprop_param_, patch_infos_[p].corners_[2]) = Param(1, 1);
		mesh().property(vprop_param_, patch_infos_[p].corners_[3]) = Param(0, 1);

		for(k = 0; k < patch_vertices[p].size(); k++){
			vh_i = patch_vertices[p][k];
			i = mesh().property(vprop_lns_idx_, vh_i);

			w_ii = 0.0;
			bu[2*i] = 0.0;
			bu[2*i+1] = 0.0;

			for(voh_it = mesh().voh_iter(vh_i); voh_it; ++voh_it){
				vh_j = mesh().to_vertex_handle(voh_it);
				j = mesh().property(vprop_lns_idx_, vh_j);

				w_ij = mesh().property(eprop_w_, mesh().edge_handle(voh_it));
				w_ii += w_ij;

				if(j >= 0){
					assert(mesh().property(vprop_side_idx_, vh_i) == p);
					adj_idx = mesh().property(vprop_side_idx_, vh_j);

					if(p != adj_idx){
						int kk = -1;
						for(kk = 0; kk < 4; kk++){
							if(patch_infos_[p].adj_patch_idx_[kk] == adj_idx){
								break;
							}
						}

						assert(kk != -1);

						Scalar w00, w01, w10, w11, b[2];
						TransFunc tf = patch_infos_[p].trans_func_[kk];

						w00 = -w_ij * tf.a;
						w01 = w_ij * tf.b;
						w10 = -w_ij * tf.b;
						w11 = -w_ij * tf.a;

						b[0] = -w_ij*tf.a * tf.t[0] + w_ij * tf.b * tf.t[1];
						b[1] = -w_ij*tf.b * tf.t[0] - w_ij * tf.a * tf.t[1];


						if(w00 != Scalar(0.0)){
							fprintf(fp_ijv, "%d %d %f\n", 1+2*i, 1+2*j, w00);
						}
						if(w01 != Scalar(0.0)){
							fprintf(fp_ijv, "%d %d %f\n", 1+2*i, 2+2*j, w01);
						}
						if(w10 != Scalar(0.0)){
							fprintf(fp_ijv, "%d %d %f\n", 2+2*i, 1+2*j, w10);
						}
						if(w11 != Scalar(0.0)){
							fprintf(fp_ijv, "%d %d %f\n", 2+2*i, 2+2*j, w11);
						}

						if(b[0] != Scalar(0.0)){
							bu[2*i] += b[0];
						}
						if(b[1] != Scalar(0.0)){
							bu[2*i+1] += b[1];
						}
					}
					else{
						fprintf(fp_ijv, "%d %d %f\n", 1+2*i, 1+2*j, -w_ij);
						fprintf(fp_ijv, "%d %d %f\n", 2+2*i, 2+2*j, -w_ij);
					}
				}
				else{
					Param uv = mesh().property(vprop_param_, vh_j);

					bu[2*i] += w_ij*uv[0];
					bu[2*i+1] += w_ij*uv[1];

					assert(vh_j == patch_infos_[p].corners_[0] || vh_j == patch_infos_[p].corners_[1] ||
						vh_j == patch_infos_[p].corners_[2] || vh_j == patch_infos_[p].corners_[3]);
				}
			}
			fprintf(fp_ijv, "%d %d %f\n", 1+2*i, 1+2*i, w_ii);
			fprintf(fp_ijv, "%d %d %f\n", 2+2*i, 2+2*i, w_ii);
		}

	}
	fclose(fp_ijv);

	taucs_logfile("stdout");


	A = taucs_ccs_read_ijv(ijv_filename, TAUCS_SYMMETRIC | TAUCS_DOUBLE);
	u = (taucs_double *) taucs_vec_create(A->n, A->flags);

	printf("linear system: %d\n", A->n);

	for (i=0; i<n*2; ++i)   { u[i] = 0; }

	taucs_ccs_write_ijv(A, "test.ijv");     // just for debugging


	/////////////////////////////////////
	//// cg solver (it is working also!)
	/////////////////////////////////////
	//printf("cg solver\n");
	//
	//for (i=0; i<n; ++i)   { printf("bu[%3d] = %lf\n", i, bu[i]); }
	//for (i=0; i<n; ++i)   { printf("bv[%3d] = %lf\n", i, bv[i]); }
	//for (i=0; i<n; ++i)   { printf("bw[%3d] = %lf\n", i, bw[i]); }
	//
	//// computing u
	//rc = taucs_conjugate_gradients(A, NULL, NULL, u, bu, 2*A->n, 1e-7);
	//if (rc == TAUCS_SUCCESS)  printf("cg solver for computing u: it's ok\n");
	//else                      printf("cg solver for computing u: it's NOT ok\n");
	//// computing v
	//rc = taucs_conjugate_gradients(A, NULL, NULL, v, bv, 2*A->n, 1e-7);
	//if (rc == TAUCS_SUCCESS)  printf("cg solver for computing v: it's ok\n");
	//else                      printf("cg solver for computing v: it's NOT ok\n");
	//// computing w
	//rc = taucs_conjugate_gradients(A, NULL, NULL, w, bw, 2*A->n, 1e-7);
	//if (rc == TAUCS_SUCCESS)  printf("cg solver for computing w: it's ok\n");
	//else                      printf("cg solver for computing w: it's NOT ok\n");


	//////////////////////////////////////////////////////////////////////////
	//// LLT factor, solve, free the factorization
	//////////////////////////////////////////////////////////////////////////
	char* options_factor[] = {"taucs.factor.LLT=true", NULL };
	void* F = NULL;
	rc = taucs_linsolve(A, &F, 0, NULL, NULL, options_factor, NULL);    // factor
	if (rc != TAUCS_SUCCESS)  printf("problem in factor\n");
	rc = taucs_linsolve(A, &F, 1, u, bu, options_factor, NULL);        // computing u
	if (rc != TAUCS_SUCCESS)  printf("problem in solving x0\n");
	//rc = taucs_linsolve(A, &F, 1, v, bv, options_factor, NULL);        // computing v
	//if (rc != TAUCS_SUCCESS)  printf("problem in solving x1\n");
	rc = taucs_linsolve(NULL, &F, 0, NULL, NULL, NULL, NULL);
	if (rc != TAUCS_SUCCESS)  printf("problem in free\n");


	// get results
	FILE *tf = fopen("param.txt", "w");
	for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)  
	{
		i = mesh().property(vprop_lns_idx_, v_it);
		if (i < 0)
			continue;

		Param   uv;

		uv[0] = (Scalar) u[2*i];
		uv[1] = (Scalar) u[2*i+1];

		mesh().property(vprop_param_, v_it) = uv;
		fprintf(tf, "%d %lf %lf\n", i, uv[0], uv[1]);
	}

	fclose(tf);

	return  0;
}

void FnSurfaceGlobalParam::swap_param_values()
{
	Mesh::VertexIter v_it;
	Param uv;
	
	for(v_it = mesh().vertices_begin(); v_it != mesh().vertices_end(); ++v_it){
		uv = mesh().property(vprop_param_, v_it);
		if(uv[0] < 0.0 || uv[0] > 1.0 || uv[1] < 0.0 || uv[1] > 1.0)
			swap_param_value(v_it);
	}
}

void FnSurfaceGlobalParam::swap_param_value(const Mesh::VertexHandle _vh)
{
	int p = mesh().property(vprop_side_idx_, _vh);

	int idx, idx2;
	TransFunc tf;
	Param new_uv, uv, new_uv2;

	uv = mesh().property(vprop_param_, _vh);

	for(idx = 0; idx < 4; idx++){
		tf = patch_infos_[p].trans_func_[idx];

		new_uv[0] = tf.a*uv[0] + tf.b*uv[1] + tf.t[0];
		new_uv[1] = -tf.b*uv[0] + tf.a*uv[1] + tf.t[1];

		if(new_uv[0] >= 0.0 && new_uv[0] <= 1.0 && new_uv[1] >= 0.0 && new_uv[1] <= 1.0){
			mesh().property(vprop_param_, _vh) = new_uv;
			mesh().property(vprop_side_idx_, _vh) = patch_infos_[p].adj_patch_idx_[idx];
			break;
		}
	}

	if(idx == 4){
		int p2;
		TransFunc tf2;
		bool flag = false;

		for(idx = 0; idx < 4; idx++){
			tf = patch_infos_[p].trans_func_[idx];
			p2 = patch_infos_[p].adj_patch_idx_[idx];

			new_uv[0] = tf.a*uv[0] + tf.b*uv[1] + tf.t[0];
			new_uv[1] = -tf.b*uv[0] + tf.a*uv[1] + tf.t[1];

			for(idx2 = 0; idx2 < 4; idx2++){
				tf2 = patch_infos_[p2].trans_func_[idx2];

				new_uv2[0] = tf2.a*new_uv[0] + tf2.b*new_uv[1] + tf2.t[0];
				new_uv2[1] = -tf2.b*new_uv[0] + tf2.a*new_uv[1] + tf2.t[1];

				if(new_uv2[0] >= 0.0 && new_uv2[0] <= 1.0 && new_uv2[1] >= 0.0 && new_uv2[1] <= 1.0){
					mesh().property(vprop_param_, _vh) = new_uv2;
					mesh().property(vprop_side_idx_, _vh) = patch_infos_[p2].adj_patch_idx_[idx2];
					flag = true;
					break;
				}
			}

			if(flag == true)
				break;
		}

		if(flag == false){
			printf("!");
		}
	}
}

void FnSurfaceGlobalParam::read_face_side_info(const std::string& _filename)
{
  std::vector<std::string>  face_group_names;
  std::size_t   i;
  FILE          *fp;

  fp = fopen(_filename.c_str(), "r");

  char    s[1024];
  int     group_idx = -1;
  int     face_idx = -1;

  memset(&s, 0, 1024);
  
  while(fp && !feof(fp) && fgets(s, 1024, fp))
  {     
    // skip comment
    if (s[0] == '#' || isspace(s[0])) continue;   

    // group
    else if (strncmp(s, "g ", 2) == 0)
    {
      char  buffer[512];

      memset(&buffer, 0, 512);
      sscanf(s, "g %s", buffer);
      
      std::string   name(buffer);
      if (name == "<noname>")
      {
        group_idx = -1;
      }
      else
      {
        group_idx = -1;
        for (i=0; i<face_group_names.size(); ++i)
        {
          if (name == face_group_names[i])
          {
            group_idx = (int) i;
            break;
          }
        }
        if (group_idx < 0)
        {
          group_idx = (int) face_group_names.size();
          face_group_names.push_back(name);
        }
      }
      //std::cout << "group " << name << " " << group_idx << std::endl;
    }
    // face
    else if (strncmp(s, "f ", 2) == 0)
    {      
      Mesh::FaceHandle    fh;

      ++face_idx;
      //std::cout << face_idx << " " << group_idx << std::endl;

      fh = mesh().face_handle(face_idx);
      mesh().property(fprop_side_idx_, fh) = group_idx;
    }

    memset(&s, 0, 200);
  }
  fclose(fp);


  // automatically find corner vertices by using group name
  int   nx_face_side_idx = -1, px_face_side_idx = -1;
  int   ny_face_side_idx = -1, py_face_side_idx = -1;
  int   nz_face_side_idx = -1, pz_face_side_idx = -1;
  for (i=0; i<face_group_names.size(); ++i)
  {
    std::string& group_name = face_group_names[i];

    if (group_name == "-x")         nx_face_side_idx = (int) i;
    else if (group_name == "+x")    px_face_side_idx = (int) i;
    else if (group_name == "-y")    ny_face_side_idx = (int) i;
    else if (group_name == "+y")    py_face_side_idx = (int) i;
    else if (group_name == "-z")    nz_face_side_idx = (int) i;
    else if (group_name == "+z")    pz_face_side_idx = (int) i;
  }

  if (nx_face_side_idx < 0 || px_face_side_idx < 0 ||
      ny_face_side_idx < 0 || py_face_side_idx < 0 ||
      nz_face_side_idx < 0 || pz_face_side_idx < 0)
    return;

  std::cout << "nx_face_side_idx: " << nx_face_side_idx << std::endl;
  std::cout << "px_face_side_idx: " << px_face_side_idx << std::endl;
  std::cout << "ny_face_side_idx: " << ny_face_side_idx << std::endl;
  std::cout << "py_face_side_idx: " << py_face_side_idx << std::endl;
  std::cout << "nz_face_side_idx: " << nz_face_side_idx << std::endl;
  std::cout << "pz_face_side_idx: " << pz_face_side_idx << std::endl;

  FnSurfaceGlobalParam::nx_face_side_idx = nx_face_side_idx;
  FnSurfaceGlobalParam::px_face_side_idx = px_face_side_idx;
  FnSurfaceGlobalParam::ny_face_side_idx = ny_face_side_idx;
  FnSurfaceGlobalParam::py_face_side_idx = py_face_side_idx;
  FnSurfaceGlobalParam::nz_face_side_idx = nz_face_side_idx;
  FnSurfaceGlobalParam::pz_face_side_idx = pz_face_side_idx;


  Mesh::VertexIter      v_it;
  Mesh::VertexFaceIter  vf_it;
  int                   face_side_idx;
  //Param                 uvw;

//  std::cout << "corner vertex infos" << std::endl;
//  for (v_it=mesh().vertices_begin(); v_it!=mesh().vertices_end(); ++v_it)
//  {
//    std::set<int>   adj_face_side_idx;
//    for (vf_it=mesh().vf_iter(v_it); vf_it; ++vf_it)
//    {
//      face_side_idx = mesh().property(fprop_side_idx_, vf_it);
//      if (face_side_idx < 0)
//        continue;
//      adj_face_side_idx.insert(face_side_idx);
//    }
//
//    if (adj_face_side_idx.size() != 3)
//      continue;
//
//   /* if (adj_face_side_idx.find(nx_face_side_idx) != adj_face_side_idx.end() &&          // -x
//        adj_face_side_idx.find(ny_face_side_idx) != adj_face_side_idx.end() &&          // -y
//        adj_face_side_idx.find(nz_face_side_idx) != adj_face_side_idx.end())            // -z
//      uvw = Param(0,0,0);
//    else if (adj_face_side_idx.find(px_face_side_idx) != adj_face_side_idx.end() &&     // +x
//             adj_face_side_idx.find(ny_face_side_idx) != adj_face_side_idx.end() &&     // -y
//             adj_face_side_idx.find(nz_face_side_idx) != adj_face_side_idx.end())       // -z
//      uvw = Param(1,0,0);
//    else if (adj_face_side_idx.find(px_face_side_idx) != adj_face_side_idx.end() &&     // +x
//             adj_face_side_idx.find(py_face_side_idx) != adj_face_side_idx.end() &&     // +y
//             adj_face_side_idx.find(nz_face_side_idx) != adj_face_side_idx.end())       // -z
//      uvw = Param(1,1,0);
//    else if (adj_face_side_idx.find(nx_face_side_idx) != adj_face_side_idx.end() &&     // -x
//             adj_face_side_idx.find(py_face_side_idx) != adj_face_side_idx.end() &&     // +y
//             adj_face_side_idx.find(nz_face_side_idx) != adj_face_side_idx.end())       // -z
//      uvw = Param(0,1,0);
//    else if (adj_face_side_idx.find(nx_face_side_idx) != adj_face_side_idx.end() &&     // -x
//             adj_face_side_idx.find(ny_face_side_idx) != adj_face_side_idx.end() &&     // -y
//             adj_face_side_idx.find(pz_face_side_idx) != adj_face_side_idx.end())       // +z
//      uvw = Param(0,0,1);
//    else if (adj_face_side_idx.find(px_face_side_idx) != adj_face_side_idx.end() &&     // +x
//             adj_face_side_idx.find(ny_face_side_idx) != adj_face_side_idx.end() &&     // -y
//             adj_face_side_idx.find(pz_face_side_idx) != adj_face_side_idx.end())       // +z
//      uvw = Param(1,0,1);
//    else if (adj_face_side_idx.find(px_face_side_idx) != adj_face_side_idx.end() &&     // +x
//             adj_face_side_idx.find(py_face_side_idx) != adj_face_side_idx.end() &&     // +y
//             adj_face_side_idx.find(pz_face_side_idx) != adj_face_side_idx.end())       // +z
//      uvw = Param(1,1,1);
//    else if (adj_face_side_idx.find(nx_face_side_idx) != adj_face_side_idx.end() &&     // -x
//             adj_face_side_idx.find(py_face_side_idx) != adj_face_side_idx.end() &&     // +y
//             adj_face_side_idx.find(pz_face_side_idx) != adj_face_side_idx.end())       // +z
//      uvw = Param(0,1,1);
//*/
//    mesh().property(vprop_type_, v_it) = kCornerVertex;
//  //  mesh().property(vprop_param_, v_it) = uvw;
//
//    std::cout << v_it << std::endl;
//	std::cout << "   is_boudary: " << (mesh().is_boundary(v_it) ? "true" : "false") << std::endl;
//    std::cout << "   pos: " << mesh().point(v_it) << std::endl;
//    //std::cout << "   uvw: " << mesh().property(vprop_param_, v_it) << std::endl;
//  }
}

